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ABSTRACT 


Ine report summarizes the research undertaken at the University of 
California Davis from 1977-1980 aimed at improving present 
capabilities for computer simulation of turbulent recirculating 
flows. Attention has been limited to two-dimensional flows and 
principally to statistically stationary motion. 

The work has been of two types: research aimed at improving the 

turbulence model and work on the development of the numerical 
solution procedure. The research on turbulence modelling has 
explored separately the treatment of the near-wall sublayer and 
the exterior (fully turbulent region) working within the framework 
of turbulence closures requiring the solution of transport 
equations for the turbulence energy and its dissipation rate. The 
work on the numerical procedure, which has been based on the 
Gosman-Pun program TEACH, has addressed the problems of 
incorporating the turbulence model as well as such matters as the 
extension to time-dependent flows, the incorporation of a 3rd- 
order approximation of convective transport and the treatment of 
non-orthogonal boundaries. 

Most of the work has already been documented in the open 
literature either as journal articles or as computer program 
guides. The present report does not attempt to duplicate these 
documents but instead provides a summarizing review of the work 
which may serve both to capture the flavor of the project as a 
whole and to provide a guide to the literature for those intending 
a deeper examination. 



1. Introduction 


Tho work summarized in this report covers the period August 1977- 
January 31st 1981. The research has been aimed at developing a 
reliable, tested method for calculating turbulent tlows involving 
separated-f low regions. Although the study has been confined to 
two-dimensional situations several different questions have been 
addressed. Broadly these fall into one of two categories relating 
either to how the turbulent mixing should be approximated (or 
“modelled") or to how the posed set of differential equations 
could be accurately solved. Section 2 of this report considers 
the latter question while Section 3 deals with the former. No 
attempt is made to duplicate already documented workj where the 
research results are published only an outline summary is 
provided. Passages of the report that contain new (l.e. hitherto 
unpublished) results or which provide a new perspective on the 
published information are presented more fully. 

J. Numerical Procedure for Elliptic Flow 


1 The Starting Point and a Summary of the Work 

The computations of recirculating flows have been founded on the 
Gosman-Pun two-dimensional, steady, elliptic solver TEACH. In its 
basic form it solves, by finite volume discretization, the 
Reynolds equations for plane two-dimensional or axisymmetric flow 
employing primitive variables, using an upwind approximation of 
convective transport and the SIMPLE algorithm (Patankar [1]). The 
turbulence model built into the code is the high-Reynolds-number 
form of the Jones-Launder [2] )c * c Boussinesq viscosity model 
matched to wall functions broadly (but not entirely) in accord 
with the proposals of Launder and Spalding [3]. 
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In the work undertaken for NASA the following corrections, 
adaptations and exten5;ions to the basic TEACH code have been 
introduced. 


(i) 


The original TEACH Code was inconsistent in the 
assumptions made about the relative location of the 
velocity and pressure nodes. The velocity nodes are 
defined as lying midway between pressure nodes. The 
program also presumes, conversely, that the pressure 
nodes are situated mid-way between velocity nodes. 

For a uniform mesh the converse statement is obviously 
true but not otherwise. 



pressure nodes 

x-directed velocity 
nodes 


The program was adapted to remove the presumption. It 
turned out that this modification did not significantly 
modify the results obtained, at any rate for the fairly 
modest expansion ratios (- 1.2:1) employed. 

(ii) Although TEACH was set up for the backward-step 

geometry we found a number of inconsistencies in the 
treatment of corner control volumes where, in certain 
cases, flow can occur through only part of a cell. 


Note: flow occurs through only a 

portion of lower face of 
control volume and wall stress 
is likewise applied to only 
a portion of this face. 
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This modification h.ul only a small influence on the 
computed flow field f r the backward step though Durst 
and Rastogi [4] found important differences for flow over 
an isolated rectangula’' sectioned block resting on a plane 
surface. (In the latter case there are strong y-direction 
velocities just in front of the rib and the x-direction 
velocity is also large) . 

(iii) Adaptations to incorporate the modified set of wall 
functions due to Chieng and Launder [13]. 

(iv) Adaptations to include an algebraic stress model (ASM) 
as an alternative to the k •• c BVM. 

(V) Inclusion of time-dependent terms. 

(vi) Adaptation to allow computation of confined backward- 
facing steps with non parallel walls. 

The ensemble of changes noted above are incorporated in a program 
VAST-STEP (Viscous and Algebraic-Stress Turbulent Simulations with 
a Teach Elliptic program) . 

We summarize briefly in the following sections the work entailed 
in (iii) - (vi). More comprehensive accounts are given by Sindir 
(51, 16]. 

Before turning to these topics we mention tvo items which* though 
explored during our research* were finally not included in the 
computer program for elliptic flows. As the program of work 
neared completion* papers began to appear which made it clear 
that* in a number of recirculating flows* the use of upwind 
differencing of convective terms led to unacceptably large 
numerical errors. The use of quadratic upstream differencing of 
Leonard [7] * formally of 3rd order accuracy* was found to give 

uniformly better results. An initial attempt was therefore 
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made to include this form following the recommendations of Han ct 
al . (8]. It was found however that for turbulent flow, stable 

results could not be obtained. Now, Han et al. [8] had shown for 
the abrupt pipe expansion (a very similar flow to the backstop) 
that upwind and quadratic-upwind differencing produced nearly the 
same results (the reason being, apparently, that these are "thin" 
recirculating flows where, apart from a small patch around the 
reattachment point, the mean streamlines are oriented at only a 
small angle with the mesh thus keeping false diffusion to 
unimportant levels*. Accordingly, since there seemed no reason to 
doubt that the numerical results obtained in our study with a 42 x 
42 mesh were sensibly grid independent (6] the original upwind 
approximation was retained.** 

In the final year of the project a "multiscale" version of the 
backstop code was created by Professor Hanjalic. The 
incorporation of this model led to essentially the same computed 
behaviour as with the single-scale scheme (due to the fact that 
nearly all the turbulent energy finds itself in the "production" 
range - a problem developed and explained in Section 3) . For this 
reason therefore no formal 'guide' to this version has been 
developed . 


*The false-diffusion coefficient is generallv estimated as 
proportional to the sine of twice the anqle*nade by the flow 
with respect to the grid line. 

••Subsequently, the winter's research grouo has entirely abandoned 
upwind dif .erencina. of the various alternatives tested (power- 
law differencinq 19], skew unwind differencinq [10], locallv 
analytic differencing scheme [1 1 ] ) the quadratic differencinq 
appears to give, overall, the best results [12]. 



2.2 Incorporation of New Wall Functions 


The term "wall function" refers to the special formulae applied in 
the 'finite-volume' equations for the wall-adjacent control 
volumes for the turbulence variables and the velocity component 
parallel to the wall. The wall functions attempt to include, in 
rather simple form, the physics of the viscosity-affected sublayer 
and near-wall region where properties vary so rapidly with 
distance. Strictly wall functions ought perhaps to be discussed 
under section 3. They are presented here, however, because the 

underlying physical model of the wall functions developed by 
Chieng and Launder [13] is the same as the standard version [3]; 
it is in the way the physics are translated to formulae for the 
near-wall cells that differences emerge. 

The principal differences occur in the balance equation for 
kinetic energy. The kinetic energy level near the wall is 

dominated by generation and destruction processes which appear as 
sources and sinks in the conservation equation . In local 

equilibrium they are the only terms of any significance and, 
outside the viscous layer, they both vary as the reciprocal of 
distance from the wall. Since a finite volume method obtains the 
•difference' equations by integration of the conservation equation 
over the control volume it follows that cell-averaged values of 
the production and dissipation are needed. In the original TEACH 
wall function a cell-averaged value of production is obtained but 
a point (i.e. the nodal) value of energy dissipation rate is used. 
This inconsistent treatment is clearly undesirable. A further 
serious (though, as it turns out, partly compensating) error is 
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that the original wall functions did not distinguish the different 
physical significance of the total shear stress, t» and the 
turbulent shear stress - puv. In the viscous layer the former 
stays finite but the latter goes to zero. Now, apart from minor 
terms, the turbulence energy generation rate for a point near the 
wall is - puv3U/3y which clearlv vanishes when uv becomes zero. 
Because they used t in place of puv the original wall functions 
greatly over estimated the turbulence energy generation in the 
near-wall cell (typically by a factor of 3 or 4) . 

The Chieng-Launder (13] wall functions incorporate a consistent 
cell averaged treatment for both generation and dissipation 
including, for the latter, an exact expression for the rate of 
energy dissipation in the viscous sublayer. The other difference 
between the original TEACH wall functions and those incorporated 
in VAST STEP is that the charccteristic turbulent velocity scale 
for the near wall region is taken as that at the outer edge of the 
viscous sublayer rather than that at the node of the near-wall 
cell. The latter choice was physically undesirable because it 
meant that the characteristic velocity scale depended (greatly, in 
the region of reattachment) on the size of the near-wall cell. In 
the Chieng-Launder scheme the kinetic energy at the sublayer edge 
was found by linear extrapolation of the values at the two nodes 
closest to the wall. It must be admitted that this modification, 
while logical, did not, from a practical point of view, reduce the 
sensitivity of the near-wall flow properties to the size of the 
near-wall cell; indeed it rendered the results somewhat more 
sensitive. The reasons for this are numerous though seem to be 
principally due to the assumed logarithmic law for velocity from 
which the wall shear stress is obtained. Research continues on 
this problem at UMIST and interim results are documented in [14]. 



I 


7 

uJ ' j 

2.3 Incorporation of Algebraic Stress Model 

The introduction of the algebraic stress model due to Gibson and 
Launder [15] (omitting buoyant terms) in place of the Boussinesq 
viscosity stress strain law initially caused divergence of the 
solution procedure. Several approaches were tried before a 
convergent treatment was evolved. The initial practice of putting 
all terms relating to the Reynolds stresses into one of two source 
arrays proved spectacularly unsuccessful. The next step was to 
evaluate the Reynolds stresses on a staggered array of points as 
shown in the diagram (for simplicity, for the case of a uniform 
mesh) . 



The staggering ensured that the shear stresses were evaluated at 
exactly the positions they were needed when formulating a cell 
momentum balance for both the U and V cells. This arrangement, 
however, still did not produce a convergent field when used in 

connection with the source terms. Finally the following 
successful iteration sequence was developed. Firstly the four 
dimensionless stress ratios, u^A, uv/k etc. were obtained 
simultaneously (and at the same points) by solving the 4x4 

matrix. These ratios were used as the ba^is for finding the new 

stresses at the next cycle of iteration. In obtaining these 

values as much as possible of the stress-strain connection was 
absorbed into an effective viscosity. For example, in a free 

shear flow the form of the ASM adopted gives: 

^Cj- 1 + P/el € [k 3y k 3xJ 


-uv 
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Thus, in analysing the x-momentum equation, we treat 

I 1 - <^2 \ 

^ci - 1 + P/c) J 

as an effective viscosity and the remainder of the stress strain 
relation, involving 3V/3x, is absorbed as a source. By proceeding 
in this way unerringly convergent iteration was achieved. 


2. 4 Time-Dependent Explorations 

After an initial phase of exploration with the elliptic program we 
asked ourselves whether the turbulent flow downstream of a 
backstop should really be treated as steady. There seemed to us 
two possible situations that would necessitate a time dependent 
approach: first, the obvious case where the flow was periodic and 

second the case where the final state was steady but where two or 
more solutions existed for the given boundary conditions and 
where, therefore, we had to march through a true time transient to 
discover which solution 'Nature* chose. 

Our concern was stimulated by laser Doppler studies of laminar 
flow in a double-step expansion (16], which showed both 
asymmetries and ti.ne dependencies above certain limiting Reynolds 
numbers. There was also felt to be a curious behaviour in the 
heat transfer data downstream of an abrupt expansion reported by 
Zemanick and Dougall [17]. For a diameter ratio of 0.43:1 the 
maximum surface heat transfer rate (that occurs close to the 
position of reattachment) appeared to occur two diameters further 
downstream at a Reynolds number of 40000 than at 20000. Such a 
change suggested some definite alteration in the flow p'jttern; 
the calculated distribution of Nusselt number, however, was 
essentially the same at the two values of R^. 


After discussions with the sponsor it was agreed that some exp- 



9 


loration of the uniqueness and/or periodic nature of the backstop 
flow be investigated for up to six months. The code was initially 
transformed to an implicit time-dependent solver by sim- 
replacing by ^ and inserting into the appro'''.'. ace 

source arrays; all other terms in the equation were evaluated r.t 
the 'new' time level. Following the experience of a group at G.E, 
Schenectady (Dr M Lubert, private communication) that this 
implicit treatment entirely suppressed the formation of a Karman 
vortex street behind a square sectioned cylinder, the code was 
adapted so that the terms other than the time-dependent one could 
be evaluated at any fraction « between the new time state (« = 1 ; 
implicit) and the old (» - 0; explicit). The adaptations required 
to the code are entirely straightforward but substantial and are 
documented in detail by Sindir [5], 16). In practice the value 
«=0.5was adopted for all the subsequent calculations. 

There followed several months of effort searching for periodic and 
or asymmetric flow behaviour in the double backstep flow. The 
approach followed was to start from rest and gradually raise to 
its asymptotic level the inlet velocity (several channel heights 
upstream of the backstep) . The resultant flow certainly showed a 
degree of oscillation but the amplitude of these oscillations 
diminished with time. A variant was to apply a slightly 
asymmetric velocity at inlet thereby producing a correspondingly 
non-symmetric flow field in the interior. When this had reached a 
stationary state the asymmetry in the inlet profile was removed 
and the computation continued to see whether a symmetric flow was 
re-established. On one occasion a spectacularly asymmetric flow 
had been established which remained after the symmetric boundary 
conditions were restored. Because of the relatively coarse mesh 
used however the flow evidently suffered from false diffusion. 
Subsequent attempts with a finer mesh did not succeed in 
reproducing such a flow. 
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After the budgeted period had been spent it was felt that the 
questions of non-uniqueness and periodicity were ones that could 
easily absorb years of work and a very large computer budget. We 
could afford neither and so we decided to proceed on the 
assumption that a steady flow model was a reasonable one for 
examining the backward facing step. (Both time-dependent and 
steady flow versions of VAST STEP are listed in [6], hc'ever.) 

As a footnote to this part of our study we note that three groups 
contributed computations of the turbulent backstop to the 1981 
Stanford Conference using time-dependent procedures. Ha Minh et 
al . [18] and Mansour and Morel [19] using the k -• e Boussinesq 
turbulence model both obtained steady flow results. Mellor and 
Celenligil [20] using a Reynolds stress closure, however, found a 
periodic behaviour. Clearly the question is still open. 

2.5 Adaptation to allow computation of non-rectllinear walls 

We desired to extend the capabilities of the numerical procedure 
to allow the computation of separated flow in weakly converging or 
diverging passages. Three possibilities were considered for 
achieving this: a) an orthogonal curvi-linear coordinate system; 

b) a non-orthogonal mesh in which one set of coordinate lines was 
straight parallel lines at right angles to the main flow; c) a 
Cartesian grid where special balance equations were developed for 
the edge control volumes which, along at least one of the walls 
were of complex form as indicated in the sketch. 


I 
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Having regard for the fact that an algebraically cumbersome 
turbulence model was to be used (the ASM) it quickly seemed that 
the last option v;as the best for that retained the relative 
simplicity of the Cartesian mesh over the interior flow. Howe^'er 
we met no success in obtaining convergent behaviou*' even for roof 
angles as small as 4 degrees. Finally, follov/ing a suggestion by 
Patankar (9) the sloping roof w< ' treated as simply a succession 
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This simplified approach certainly converged for both ASM and BVM 
treatments and it is this which has been built nto VAST STEP. It 
appears to be satisfactory provided that the angles of exparsion 
are small (no more than about 8 degrees) and provided that the 
flow region of interest is not along the "stepped" wall itself. 
These conditions were met in the Stanford Test Cases that were 
tried. The view that the approach was "satisfactory" is based on 
the fact that the computations [6], [21] shew a similar 
disagreement with the data as were shewn for the case of parallel 
w*i Is, 

The limitations of the method are evident however and there is 
clearly a need, in the near future, to embed an ASM into an 
arbitrary orthogonal coordinate system. 
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3. Turbulence Modelling Research 


3.1 The Overall Plan t 

In the first year of work attention was confined to thin shear 
flows. Moreover only the fully turbulent regiv^n was considered 
for in the boundary layers in question the near-wall viscosity- 
affected region was well described by the "universal" semi- 
logarithmic velocity profile. The work was initially aimed at 
placing on a firmer footing than hitherto the "multi-scale" 
modelling approach which had been developed by Launder & Schiestel 
[22, 23]. Section 3.2 outlines the principal feati’.es of this 

work. It was discovered, however, that most of the benefits that 
we felt accrued from the multi-scale treatment could in fact -.Iso 
be obtained with a single i cale treatment provided appropriate 
modifications were introduced to the dissipation rate equation. 
In the second year attention shifted to turbulent separated flows. 
The handling of the near-wall region under conditions far removed 
from local equilibrium via wall functions has already been 
discussed in Section 2.2, Likewise, work on incorporating an 
algebraic stress closure into the recirculating flow procedure 
posed mainly numerical problems and was discussed in Section 2.3. 
The third year's work combined the exploration of new closure 
ideas and their evaluation in both recirculating am thin shear 
flows. Several faults of the turbulence model were identified and 
attempts made to remove them; these are discussed in Section 3. 
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3 . 2 Tho Multi-Pcnlo Turbulence Closure - and its Single scale 

derivative 


The multi-scale approach was evolved to cope with the fact that 
the local rate of energy dissipation, c , is not as directly 
connected with the mean flow anu energetic motions as any of the 
conventional rate equations for f. suggested. The energy 
containing part of the spectrum was thus divided into two portions 
as indicated in figure 3.1: a production region and a transfer 

region. This move naturally brought into focus the quantity Cp, 
the transfer rate of energy from the former region to the latter. 
The multi-scale model of Launder and Schiestel [22, 23] provided 

independent transport equations for finding both Cp and 
(quantities that would be equal if spectral equilibrium 
prevailed) . While ‘ p is associated with large-scale interactions 
and is directly responsive to the mean strain field, e is held to 
feel these influences only indirectly through their effect on the 
energy transfer rate, e^. These ideas are broadly in accord with 
the accepted view of the dynamics of the turbulent energy 
spectrum. 

Work undertaken by Dr Hanjalic during his period of support as a 
visiting research engineer was directed at improving the 
generality of this model. While the framework of the model 
remained intact, the detailed form underwent considerable change 
in order to provide, so far as possible, an internally consistent 
physical picture. Moreover the range of flows was extended from 
the three equilibrium shear flows considered by Launder and 
Schiestel to include boundary layers and distorted grid 
turbulence. The out-come of this study is •'^ported in Hanjalic, 
Launder and Schiestel [24]. 


1 
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Towards the end of the development work on the multi-scale 
approach it was found (Hanjalic and Launder [25]) that practically 
as good results were obtainable with the single scale scheme if 
corresponding modifications to those devised for the multi-scale 
closure were introduced to the (single) equation for the energy 
transfer/dissipation rate. Only for the case of strongly acceler- 
ated grid turbulence was there a decisive advantage from using the 
multi-scale model and this example is far removed in kind from 
most practically interesting shear flows. 

The failure to achieve significantly better results with the 
multi-scale model is believed to be due to the fact that, as 
presently constructed, the energy in the "transfer" range amounts 
to only about 10% of the total kinetic energy for a typically 
shear flow. There is thus only marginal advantage to be gained 
from recognizing the difference between the rate that energy is 
fed to this transfer range and the rate that it leaves it through 
viscous dissipation. The best way forward appears to be to re- 
formulate the model placing a greater proportion of the total 
energy into the transfer range. This would necessitate that the 
assumption of isotropic turbulence in the transfer range be 
discarded, however; instead, an algebraic stress model for the 
transfer-range stresses would be appropriate. 

3. 3 Final Explorations in Turbulence Modelling 


By late 1979 we had thus concluded that provided we introduced a 
new source term proportional to the square of the mean vorticity 
(24, 25] we could achieve essentially as good agreement in thin 
shear flows with a single-scale model as with the multi-scale 
scheme. The new term had a vital effect for it made the flow more 
sensitive to normal-strain production than shear strain. The 
modification allowed much improved prediction of the axisv'mmetric 
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jet and the turbulent boundary layer in adverse pressure 
gradients. In each case, due to the flow deceleration, the levels 
of viscosity were reduced producing, for the former, lower rates 
of spread and, for the latter, a diminution in wall friction [251. 
It was therefore decided to make this single-scale model the basis 
for the recirculating flow studies. 

When, however, the small-modifications to the e equation were 
introduced into the backstep code (which, at that time, contained 
only the BVM treatment) it was found, against all expectations, 
that the computed reattachment length was substantially shortened, 
indicating a pronounced increase in turbulent viscosity. This was 
a quite unacceptable result, for the standard k- e model already 
gave significantly too short reattachment lengths. The cause was 
traced to the fact that, due to the invariant form of the new 
source term it had the unfortunate property of strongly modifying 
e levels in flows with streamline curvature. The sign of the 
effect was the opposite of what is found to occur in practice, 
e.g. when the fluid's angular momentum increased with radius 
(stable stratification) the model strongly enhanced viscosity 
levels. 

We finally concluded that streamline curvature was such an abiding 
feature of turbulent shear flows that it would be intolerable to 
retain the proposed form of dissipation equation. Instead we 
decided ’ o include a term of the same type but with the sign 
reveri.,d. In fact for simplicity the form adopted with both the 
ASM and BVM was 



Source 
of e 


( 3 . 1 ) 
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For the case of a simple shear flow the source takes the standard 
BVM form. Computations were reported using eq (3.1) for the 
various backward-step cases selected for the Stanford Conference 
[21]. The ASM computations used in conjunction with the above e 
source gave the beot agreement with experiment in the separated- 
flow region of any of the contributed results for this case. 


A further fault of the original e equation was diagnosed in the 
course of our work. In near-wall turbulence it is found that the 

3/2 

lengthscale, k /e, produced by that equation became too large in 

the vicinity of the wall as the flow approaches separation* It is 

believed to be this property which is mainly responsible for the 
failure to predict separation in various supersonic flows [26] and 
for the excessive levels of heat transfer predicted downstream of 
an abrupt pipe expansion [17] when the low-Reynolds-number version 
of the BVM is used [12]. Although this property of the e equation 
had been known for some years [3] it was only with the detailed 

measurements of East and Sawyer [27] that it became fully apparent 

that this was an undesired characteristic. Their experiments of 
equilibrium turbulent boundary layers in adverse pressure 
gradients indicated that the near-wall values of k ^/e were, to a 
close approximation, only a function of distance from the wall and 
thus independent of the turbulence energy generation to 
dissipation rate. We therefore attempted to find a further source 
term for the e equation that would counteract the tendency for the 
predicted length scale to become too large. Various source terms 
were explored involving typically spatial gradients of the kinetic 
energy. While we succeeded in devising fairly successful forms 
for boundary layer flows these uniformly led to disastrous predic- 
tions in free shear flows. An alternative approach we considered 
was to let the turbulent Prandtl number for e diffusion become a 
function of the turbulence energy generation to dissipation rate. 
An analytical study for the case where P/e was zero led to the 
idea that one might adopt the form 

2 “ 1 


1.3 


=£2 - ‘’^='=£1 


( 3 . 2 ) 
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This version did indeed load to a small overall improvement in the 
backstop computations without apparently harming the prediction of 
free flows. However the introduction of P/ e made it so difficult 
to achieve converged numerical results that we finally did not 
include the modification. 


Closing Remarks 


Despite the considerable effort expended in turbulence model 
development over the course of the research grant it must be 
acknowledged that no fully satisfactory form of length-scale 
determining equation was developed. We found no invariant way of 
removing the tendency for the near-wall length to grow too large 
without giving the equation the wrong sensitivity to streamline 
curvature. At present the writer does not feel that a satis- 
factory solution cannot be devised but clearly more thought and 
careful testing is needed. The current practice of his group is 
merely to overwrite the near-wall length scale at the value found 
under local-equilibrium conditions (i.e. P * e ) should the 
dissipation equation return length scales larger than that. We 
would be interested to learn the experience of others with such a 
modification . 


The multi-scale approach/ while of definite advantage in certain 
rapidy distorted flows/ does not/ as currently constructed/ have 
any sufficiently marked practical benefits to justify its use in 
shear flows (entailing/ as it does, the solution of two extra 
transport equations) . The potential usefulness of the approach is 
widely acknowledged but a formulation needs to be developed that 
gives a more even distribution of energy between the two regions. 
At present several groups are trying to understand the physical 
processes of energy transfer across the spectrum more quantita- 
tively as a preliminary to further model development. 
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